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As computers increase in speed and memory, scien- 
tists are inevitably led to simulate more complex systems 
over larger time and length scales. Although a simple, 
straightforward algorithm is often the most efficient for 
small system sizes, especially when the time needed to 
implement the algorithm is included, the scaling of time 
and memory with system size becomes crucial for larger 
simulations. 

In our studies of hysteresis and avalanches in a sim- 
ple model of magnetism (the random-field Ising model at 
zero temperature), we often have found it necessary to 
do very large simulations. Previous simulations were lim- 
ited to relatively small systems (up to 900 2 and 128 3 jj], 
see however [0). In our simulations we have found that 
larger systems (up to a billion spins) are crucial to ex- 
tracting accurate values of the universal critical expo- 
nents and understanding important qualitative features 
of the physics. 

We have developed two efficient and relatively 
straightforward algorithms which allow us to simulate 
these large systems. The first algorithm uses sorted 
lists and scales as 0(N log N), and asymptotically uses 
N x (sizeof(double)+sizeof(int)) bytes of memory, where 
N is the number of spins. The second algorithm, which 
does not generate the random fields, also scales in time 
as O(NlogN), but asymptotically needs only one bit of 
storage per spin, about 96 times less than the first algo- 
rithm. Using the latter algorithm, simulations of a billion 
spins can be run on a workstation with 128MB of RAM 
in a few hours. 

In this column we discuss algorithms for simulating 
the zero-temperature random-field Ising model, which is 
defined by the energy function 



Si , 



(1) 



<hj> 



where the spins Si = ±1 sit on a D-dimensional hypercu- 
bic lattice with periodic boundary conditions. The spins 
interact ferromagnetically with their z nearest neighbors 
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with strength J, and experience a uniform external field 
H(t) and a random local field hi. We choose units such 
that J = 1 . The random field hi is distributed according 
to the Gaussian distribution p{h) of width R: 
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The external field H(t) is increased arbitrarily slowly 
from — oo to oo. 

The dynamics of our model includes no thermal fluc- 
tuations: each spin flips deterministically when it can 
gain energy by doing so. That is, it flips when its local 
field 



hf = J^2s j + h i + H 
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changes sign. This change can occur in two ways: a 
spin can be triggered when one of its neighbors flips (by 
participating in an avalanche), or a spin can be triggered 
because of an increase in the external field H (t) (starting 
a new avalanche) . 

The zero-temperature random-field Ising model was 
introduced by Robbins and Ji H to study fluid inva- 
sion in porous media and front propagation in disor- 
dered systems. We have used the same model § m 
a different way |6) to model noise in hysteresis loops 
in disordered materials. In particular, we wish to un- 
derstand Barkhausen noise in magnetic materials with 
quenched disorder Q. It has been found experimentally 
that when an external field is gradually applied, many 
materials magnetize not continuously, but in a noisy way, 
with jumps (avalanches) of all sizes. (The noise can be 
heard by wrapping the magnetizing material in a coil of 
wire and amplifying the signal into a speaker. The sig- 
nal makes a crackling noise when a permanent magnet is 
brought close, quite similar to the crackling noises heard 
in fires, crisped rice cereals, and crumpled paper 0.) In 
the steepest part of the hysteresis loop, these avalanches 
are found to have a power-law distribution of sizes with 
an exponent r w f .5. Power laws are also found in the 
distribution of avalanche times with an exponent a ~ 2 
and in the power spectrum. 

The zero-temperature random-field Ising model is 
interesting because, as in the disordered magnetic ma- 
terials it attempts to model, the avalanches can have a 
broad range of sizes. If all the avalanches were small, 
understanding them would be straightforward and not 
very interesting. Indeed, at large disorder R, the chance 
that a spin which has just flipped will trigger one of its 
z neighbors scales roughly as zJ/R. If this quantity is 
smaller than unity (large disorder) , all avalanches will be 
small: the noise will be a series of small pops all of about 
the same size. This behavior is uninteresting not because 
it is simple, but because the behavior is strongly depen- 
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dent on the details of the model at short distances, where 
the model is at best a caricature of a real material. 

It also is easy to understand the system in the small 
disorder regime zJ/R >> 1, where almost all the spins 
flip over in one infinite avalanche. There are many prob- 
lems (for example, fracture and first-order phase transi- 
tions) where a single nucleation event leads to the release 
of the stored energy in a single catastrophic event. 




FIG. 1. A three-dimensional view from one side of a single 
avalanche in a 200 x 200 x 200 system at R = 2.3 (within 6% 
of the critical disorder R c ). The avalanche contains 282,785 
spins. The time when each spin flipped is shown by its color. 
The avalanche generally grew from left to right. Note that it 
has many branches and holes; the large avalanches in three 
dimensions probably have a fractal dimension a little less than 
three. Also notice that on the right hand side, there are sev- 
eral dark red spots poking through in the middle of the light 
green area. The green area stopped growing, but other parts 
of the avalanche later filled in the holes. 

We focus on the crossover between these two limit- 
ing cases, where the system exhibits crackling noise with 
avalanches of all sizes. For a particular value of the disor- 
der R = R c , a spin which has just flipped will on average 
flip exactly one neighbor as the external field H (i) is in- 
creased to a particular value H c . The avalanches at R c , 
H c (the critical point), are finely balanced between stop- 
ping and growing forever. They advance in fits and starts 
(see Figs. [l] and ||) and come in all sizes (Figs. ^ and |^) 
with a probability which decreases as a power law of the 
number of spins in the avalanche. At H c , the distribu- 
tion of avalanche sizes decays with an exponent of r w 1.6 
(quite close to the experimental results), and integrated 
over all H, the distribution decays with an exponent |T3| ] 
f 2. Below the critical disorder R c , there will be an 
avalanche which will flip a nonzero fraction of the spins 
in the system even as the system size goes to infinity: we 
call this avalanche the infinite avalanche. There are very 
large avalanches even for disorders far above the critical 



disorder. In three dimensions, there are still two decades 
of power law scaling 50% above the critical point. How- 
ever, the convergence to the expected asymptotic power 
law is very slow (Figure |^) . This behavior means that we 
see critical scaling even if we do not fine tune R to R c , 
but we must use very large systems to get close enough 
to R c to obtain a convincing power law. In practice, we 
needed simulations of approximately a billion spins to 
understand the physics in three dimensions JlO| . Two di- 
mensions remains a challenge because the proper scaling 
is not clear even for 30, 000 2 spins pC|,H. 
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FIG. 2. A time series showing the number of spins which 
flipped in each shell of the avalanche shown in Fig. i] Q . Note 
that the avalanche is a series of bursts: near the critical point, 
the avalanche is always on the verge of halting, so it proceeds 
in fits and starts. 

It is crucial with this many spins that our algo- 
rithms be efficient both in computer time and mem- 
ory. We begin by giving the simple, but inefficient 
approach which has an execution time which scales as 
0(N 2 ). We then develop a more efficient approach using 
a sorted list which gives an execution time which scales 
as 0(N log TV), but which needs memory storage which 
scales as N x (sizeof(double)+sizeof(int)). A billion spins 
would demand 12 Gigabytes of RAM for efficient execu- 
tion, which is not usually available. Finally, we give an al- 
gorithm whose execution time also scales as 0(N log N), 
but whose memory requirements are asymptotically only 
one bit per spin. In this case 10 9 spins requires 120 MB of 
storage, which is feasible on a standard workstation [0. 
We conclude with a discussion of time and space issues for 
calculating and storing histograms and correlation func- 
tions. 
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FIG. 3. Distribution of avalanche sizes for different values 
of the disorder 7? in three dimensions. Some avalanches re- 
main large (hundreds of spins) for R a factor of two above 
the critical value R c ~ 2.16 where we expect a pure power 
law. The avalanches are enormous (millions of spins) when 
the system is still 4% away from the critical point; for this 
reason we need large systems. The inset is a scaling collapse 
of the data: the thin lines in the main figure show the scaling 
prediction for the avalanche sizes stemming from the scaling 
collapse p0| . Note that the scaling predictions already work 
well at R = 4. The pure asymptotic power law behavior is 
not yet seen at R — 2.25, when six decades of scaling are 
observed. We needed simulations of a billion spins to show 
convincingly that the power law would eventually occur Jj"o[ . 




FIG. 4. A 30,000 x 30,000 simulation with disorder 
R = 0.65, where each pixel represents a 30 x 30 square, 
and each avalanche is a different color. Note that there are 
avalanches of all sizes, with many smaller avalanches, and 
fewer large ones. 



The Brute Force Method 

The brute force method is the easiest one to imple- 
ment and is competitive for system sizes up to about 
10,000 spins. In this method, we store a spin direction 
and a random field for each site of the lattice. We can 
then proceed as an experimentalist would by measuring 
the magnetization at specific predetermined values of H. 
We start with magnetization M — —N and a large neg- 
ative field Hq and then increment to H±, check all spins 
in the lattice, and flip those spins in a positive local field. 
Then we must check the neighbors of the flipped spins 
again to sec if their local fields are now positive. This pro- 
cedure is continued until all the neighbors of flipped spins 
have been checked. We then repeat the whole procedure 
again for a new field H%, and so on. This approach gives 
the correct magnetization at the fields H n : the order in 
which spins arc flipped can be shown not to influence the 
final state fL5lj5|]. However, unless the increments in H 
are very small, several avalanches may occur in a given 
increment, and all information about single avalanches 
(such as histograms of avalanche sizes) will be distorted. 

The time for the brute force method scales as 
O(NXT), where X is the number of fields H n at which 
the magnetization is measured, and T is the average time 
needed to check the neighbors of the flipped spins mea- 
sured in units of shells of neighbors. 
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FIG. 5. Example of how a queue is used to propagate an 
avalanche. The colored spins are spins which either have 
flipped in the current avalanche or will flip in the current 
avalanche. Spin 13 triggered the avalanche, then the light grey 
spins (14,8,12) were put on the queue as the first shell. As 
they flipped, the second shell, the blue spins (15,19,7,11,17), 
were put on the queue. As the first blue spins (15,19,7) 
flipped, the dark red spins (10,20,20,18,6) were added to the 
queue as the start of the third shell. The next spin to flip is 
at the left hand side of the queue. When this spin flips, its 
neighbors will be checked, and the neighbors which are ready 
to flip will be added to the right hand side of the queue. The 
small numbers below the spins in the queue indicate which 
neighbor caused the spin to be put on the queue. Note that 
different neighbors can cause a spin (such as spin 20) to be 
put on the queue more than once. We have to be careful to 
only flip the spin once. 

A variation on this approach is to propagate one 
avalanche at a time as shown in Fig. [5]: 
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1. Find the triggering spin for the next avalanche by 
checking through the lattice for the unflipped site 
with the largest internal field h\ nt — h° s — H. 

2. Increment the external field so it is just large 
enough to flip the site, and push the spin onto a 
first-in first-out (FIFO) queue (see Fig. ||, right). 

3. Pop the top spin off the queue. 

4. If the spin has not been flipped, flip it and push all 
unflipped neighbors with positive local fields onto 
the queue. 

5. While there are spins on the queue, repeat from 
step 3. 

6. Repeat from step 1 until all spins are flipped. 

This method is standard for avalanche propagation 
problems. It also is related to the propagation of cluster 
flips in the Wolff algorithm |l(| . Using a queue instead of 
recursion has two advantages. First, recursion is slower 
and more memory intensive, because each recursive call 
must push all local variables and all registers onto the 
system stack (which usually has a pre-allocated limit). If 
we use our own queue, we need only to push the coordi- 
nate of the next spin on the queue each time, and we can 
make the queue as large as necessary. Second, in order 
to produce a natural spin-flip order, we want to flip all 
spins that are ready to flip at a given time before we flip 
the spins that they cause to flip. If we put spins that are 
ready to flip on a FIFO queue, we correctly achieve this 
order. This procedure corresponds to doing a breadth- 
first search. Recursion, which is the same as putting the 
spins on a LIFO stack, would explore all possible conse- 
quences of flipping the first neighbor it looks at before 
it considers the second neighbor. This depth-first search 
produces an unnatural spin-flip order (although the final 
state after the avalanche is unchanged (T^J|]). The dy- 
namics during the avalanche of Fig. ^ assumes one shell 
of spins flipped during each time slice, which is easy to 
determine by placing markers on a FIFO spin queue, as 
shown in Fig. |^. Each time the marker is popped off of 
the queue, a new shell is started and the marker is put 
back on the end of the queue. 

Doing the brute force algorithm one avalanche at 
a time is very inefficient except at very low disorders. 
Sweeping through the entire lattice for each avalanche 
takes O(N) time per avalanche. Because there are O(N) 
avalanches, the total running time scales as 0(N 2 ). A 
hybrid approach, finite steps in field followed by internal 
propagation of avalanches, could be quite efficient if one 
is solely interested in the magnetization at those fields. 
A brute force method is probably necessary when simu- 
lating systems with long-range interactions 



Time Efficiency: Sorted Lists 

The brute force method is very inefficient at locating 
the origin of the next avalanche, and we are immediately 
led to think of storing the several largest local fields in 
each sweep. If we take this thinking to its logical conclu- 
sion, we are led to store a list of all of the spins in the 
system, sorted according to their random fields. 
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FIG. 6. Example of how a sorted list is used to find the next 
spin in the avalanche. The colors indicate spins which have 
already flipped. The first column in the sorted list contains 
the random field, and the second column contains the number 
of the spin with that random field. The arrows to the right 
indicate the nextPossible [nf] pointers — the first spin which 
would not flip with n-f neighbors up. The spins pointed to are 
the possible starting locations for the next avalanche. Note 
that some of the pointers point to spins that have already 
nipped, meaning that these spins have more than n-f neighbors 
up. In a larger system, the unflipped spins will not all be 
contiguous in the list. 

Unfortunately, life is complicated by the fact that 
spins experience not only their local random fields, but 
also fields from their neighbors. To find the origin of the 
next avalanche, we use the following algorithm: 



f. 



2. 



3. 



Define an array nextPossible [nj], = 0, 1 . . . z, 
which points to the location in the sorted list of the 
next spin which would flip if it had rif neighbors. 
Initially, all the elements of nextPossible[rjf ] point 
to the spin with the largest local random field, hi. 

Choose from the z + 1 spins pointed to by 
nextPossible, the one with the largest internal field 

tin.tr 



if] = n-f — n\ + hi = 2ri| 



Move the pointer nextPossible [rif 
on the sorted list. 



z + hi. 

to the next spin 



4. If the spin with the largest /i^ nt [n-f] has up neigh- 
bors, then flip it. Otherwise go back to step 2. 

An example of the sorted list and the pointers from 
nextPossible is shown in Fig. [| 
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The sorting of spins can be done in time 
0(N log N). Storage with this algorithm is N x 
sizeof(int) for the sorted array (if we reduce the D- 
dimensional coordinates to one number |18|), and N x 
(sizcof(spin) + sizcof( double)) for the lattice itself. Var- 
ious other compromises between speed of execution and 
storage are possible, but all leave the running time 
0(N log N). The sorted-list algorithm is fast: the largest 
system sizes we can store on a reasonable workstation 
execute 1000 2 and 100 3 spins in a few seconds. It is the 
method of choice for these small systems or when one is 
interested in the behavior for non-monotonically increas- 
ing fields. 

Space Efficiency: One Bit per Spin 

The combination of the rapid execution of the 
sorted-list algorithm and large finite size effects led us 
to develop an algorithm optimized for memory efficiency. 
The key is to recognize that we need never generate the 
random fields! In invasion percolation pl| (and in the in- 
terface problem |3| analogous to ours) the random fields 
are generated only on sites along the boundary of the 
growing cluster. In our problem, we can take this idea 
further: for each change in a spin's local field given by 
Equation |§|, we generate only the probability that it will 
flip. Storing the random fields is unnecessary because 
the external field, the configuration of the spin's neigh- 
bors, and the knowledge that the spin has not yet flipped 
gives us all the information which we need to determine 
the probability that the spin will flip. The only quan- 
tity which we must store for each site of the lattice is 
whether the spin is up or down. Thus, we can store each 
site of the lattice as a computer bit saving large amounts 
of memory. 

For a monotonically increasing external field, the 
conditional probability that a spin flips before its nonran- 
dom local field, H nl - = H + (2ra j — z), reaches H m + AH ni 
given that it has not flipped by LT nr is 



P mp (H nT ,AH nT ) = 



A(tfnr) 



(4) 



where Pj(i? nr ) is the probability that a spin points down 



*The sorted-list algorithm can be used for non-monoton- 
ically increasing fields with only a few minor additions. When 
the external field is being lowered instead of raised, the 
avalanche propagation is the same, except spins are flipped 
when their local field becomes less than zero instead of when it 
becomes greater. The nextPossible array needs to be handled 
carefully. The next spin that would flip up if the field were in- 
creased is the last spin that would have already flipped down 
with the field decreasing. Every time the direction of change 
of the external field is reversed, all of the nextPossible[n|] 
pointers need to be adjusted by one to account for this. 



when the local field is H ni . A spin with local field i? nr will 
point down if its random field hi satisfies hi + H nl - < 0. 
This condition implies that the probability that a spin 
with rif up neighbors points down is 



Pl(n h H)=J_ 



-H nT (n-f,H) 



p(h) dh 



= \+ l -evi(-H 11T {n h H)/V2R 
= ~etfc(H nx (n h H)/V2Ry 



(5) 



(6) 



(Writing P± in terms of the erfc function removed some 
problems with rounding at large negative fields H.) 
These probabilities are illustrated graphically in Fig. [?|. 
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FIG. 7. The probability that a spin will not have flipped 
by the time its local field reaches L is the probability that the 
random field is less than — L. This probability is represented 
by the shaded area of the Gaussian. The probability that the 
spin will flip before the field reaches —(L + AL) is represented 
by the area of the darker region divided by the area of the 
shaded region. 

Finding the next avalanche is subtle when the ran- 
dom fields are not stored: changing the external field H 
introduces a probability that any unflipped spin in the 
lattice may flip. Inspired by the continuous time Monte- 
Carlo algorithm |^| , we keep track of N nr , the number 
of down spins which have n-p up neighbors. Given the 
probabilities that spins with n-\ up neighbors will flip, 
we calculate both the change in the external field ALT 
needed to flip the next spin and the probability that the 
next spin to flip has rtj up neighbors. We then randomly 
choose n f and search at random through the lattice for a 
spin with n-j- up neighbors. The time taken for the search 
is the part of the algorithm which scales worst for large 
N . If there are N ny spins left, this search will take an av- 
erage time 0(N/N n ). Summing over iV„ T and yields 
a bound of order ziVlogiV. In one of our programs, we 
use a tree structure to do this search more efficiently; 
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this complication decreases the running time by 40% for 
a 500 2 system at R = 1. 

How do we calculate AHI From Equation [|, the 
probability that a single spin with n-f neighbors up has 
not flipped in the range AH is 1 — Pfu p (n^ , H, AH) = 
Pl(ni,H+AH)/Pi(nt,H). The probability that no spin 
with n-f up neighbors has flipped in this range is 



pnonc 



P^n^H + AH) 



N„ 



(7) 



and the probability that no spin has flipped between H 
and H + AH is 



P nonc (AH) = Jl Pn°' 



(8) 



To find AH, we choose a random number r uniformly 
distributed between zero and one, and set AH so that 
P nonc (AH) = r. 

Unfortunately, we cannot solve for AH analytically, 
and we must find a numerical solution. To find this solu- 
tion efficiently, we need a good initial guess. In analogy 
with nuclear decay, if spins flip with a constant rate F, 
we expect the probability that no spins have yet flipped 
to be e~ rAH . So, for a first approximation of AH, we 
assume that the spin-flip rate T is a constant, and there- 
fore 



AHi = 



log(0 
' T(H) • 



(9) 



where is given by 

dlog(P nonc (AH = 0)) 



T(H) 



dAH 

dlog(n(n T)J ff + AJ0) 



dAH 



= 2,^t p( H) =^T(n h H) (10) 

We can make a better second guess by looking at the 
error in our first guess. If the error in our guess is Ar = 
P none (AH) — r, then we can make an improved second 
guess for AH by aiming for r — Ar: 



AH, 



log(r — Ar) 



(11) 



These two guesses can then be used as input into a root 
finding routine |2~l] ]. Note that while these guesses are 
usually very good for small \H\ and lead to quick solu- 
tions, they can be very bad for large \H\. If the guesses 
for AH are very large, it may be better to choose two 



arbitrary guesses. In our code, if AHi > 20, we use 
AHi — and AH 2 — 20 for the two guesses. 

Our algorithm for finding the next avalanche be- 
comes 

1. Choose a random number r uniformly distributed 
between zero and one. 

2. Pre-calculate the values of Pi(n^,H) using Equa- 
tion ^. These values will be used repeatedly in solv- 
ing for AH. 

3. Calculate guesses for AH using Equations || and 
0, and use them as input to a root finding routine 
to find the exact solution for AH. 

4. Increment H by AH. 

5. Calculate the array probFlip[n-f ] for use in the re- 
mainder of the avalanche, where probFlip[ri|] is the 
probability at the current field H that a spin will 
flip when its number of up neighbors changes from 
raj to rif +1 (see Equation ||). 

6. Calculate the rates for flipping spins for each at 
the current field H = H \d + AH: 

T(n T , H) = N n , p(H nl (n h H))/P l (n T , H) (12) 

and the total rate T(H) = J2n 1= o r 0*Ti H )- 

7. Choose a random number uniformly distributed be- 
tween zero and T and use it to select n-\ . 

8. Search at random in the lattice for an unflipped 
spin with n-f up neighbors p2[ . 

9. Start the avalanche at that spin. During an 
avalanche, the algorithm is essentially the same as 
the brute force algorithm: 

10. Push the first spin onto the queue. 

11. Pop the top spin off of the queue. 

12. If the spin is unflipped, flip it, find rx-p, and decrease 
N nr by one. Otherwise, skip to step 14. 

13. Look at all unflipped neighbors. For each unflipped 
neighbor, find the current number of up neighbors, 
nj; decrease N ni _i by one, and increase A^„ T by 
one. Push the spin on the queue ]23[ ] with proba- 
bility probFlip[rj| — 1], as calculated in step (5). 

14. While there are spins left in the queue, repeat from 
step 11. 

15. While there are unflipped spins, repeat from step 
1. 



6 



This algorithm is about half as fast in practice as the 
sorted-list algorithm which is faster than we expected. 
The overhead involved in solving for AH is presumably 
compensated by the time saved not shifting data in and 
out of cache. Systems of 10 9 spins take a few days of 
CPU time on a reasonable workstation; 30, 000 2 systems 
take less than 15 hours on a 266 MHz Pentium II. 

Calculating Histograms and Correlations 

Several functions are needed to characterize the crit- 
ical properties of our model. The simplest function is the 
magnetization M as a function of the external field H. 
We also calculate distributions of avalanche sizes (Fig. ||, 
top), and correlation functions. Some care must be taken 
to make sure that the calculation of these functions does 
not dramatically increase the running time or memory 
requirements of the simulation. 

When doing calculations with a billion spins, we 
cannot output any quantity which scales linearly with 
the system size. Instead of computing H(M) at each 
avalanche (about a GB of data, which rapidly would fill 
our disk), we are forced to compute H(M n ) at pre-chosen 
points. 

The characteristic feature of the critical point is the 
appearance of an infinite avalanche. The equivalent of an 
infinite avalanche in a finite system is an avalanche which 
spans the entire system in at least one dimension. To 
tell whether we are above or below the critical point, we 
need to detect these spanning avalanches. In three and 
higher dimensions, the number of spanning avalanches 
as a function of R is also interesting to study. The most 
obvious way of detecting spanning avalanches is to mark 
each row as a spin flips in it, and check at the end of the 
avalanche to see if all the rows contain flipped spins from 
the avalanche. However, this method requires 0{N 1 / D ) 
operations per avalanche. Because there are many small 
avalanches, this method is unacceptable. A preferable 
method is to keep track of the 2 x D boundaries of the 
avalanche as it grows. If a pair of boundaries meet, then 
the avalanche is a spanning avalanche. We must take 
care to treat the periodic boundary conditions properly. 

Another useful function is the avalanche size distri- 
bution D(S), defined as the number of avalanches which 
flip S spins during the simulation, divided by the total 
number of spins. Like the M(H) curve, the avalanche size 
distribution scales linearly with the system size. Thus, 
we need bins up to size N, the size of the largest possible 
avalanche. Logarithmic binning is the obvious solution, 
with bin n including all sizes < S < 6™. We have 
chosen b a from 1.01 to 1.1. Large bins are preferable for 
lower statistical noise. This choice is particularly impor- 
tant in the tail of very large avalanches, where small bins 
would contain few avalanches. However, very large bins 
will systematically alter the shape of the scaling functions 
(although they will not change the critical exponents) . It 



is important to divide the final population in each bin by 
the number of integers contained within the bin (and not 
just the bin width). Clearly we should also ignore the 
early bins which do not contain any integers. 

We calculate the correlation function G(x, R) within 
an avalanche, where G(x, R) gives the probability that 
the first spin in an avalanche will cause a spin a distance 
x away to flip in the same avalanche. At the beginning 
of each avalanche, we record the coordinates of the first 
spin in the avalanche. Then, for each subsequent spin 
in the avalanche, we calculate the distance x to the first 
spin, and add one to the appropriate bin. Logarithmic 
binning is not necessary for the correlation function, be- 
cause the size of the correlation function is proportional 
to the length of the system, not the total number of spins. 
Thus, we use a fixed bin size b c = 1. At the end of the 
simulation, each bin should be normalized by the number 
of spins which are between x — b c /2 and x + b c /2 away 
from the origin. 

The only tricky part of calculating G(x, R) comes 
from the periodic boundary conditions. If the avalanche 
crosses a boundary, two points at opposite ends of the 
avalanche can come close together. Because we do not 
calculate G(x, R) for spanning avalanches, we know that 
there will be at least one row in every dimension which is 
not touched by the avalanche. To calculate separations, 
we use the periodicity of the lattice and the continuity 
of the avalanche to shift the coordinates so they are all 
on one side of these empty rows. Because we are already 
keeping track of the boundaries of the avalanche for the 
detection of spanning avalanches, finding an empty row 
is easy. 

The running times of the three algorithms as a func- 
tion of system size are shown in Fig. |8j. The brute force 
algorithm can be useful when one cares only about M(H) 
at a few points, but is otherwise too slow for large sys- 
tems. The sorted-list algorithm is the fastest algorithm, 
but on a 128 MB machine, only system sizes of about six 
million spins can be run. The bits algorithm is almost as 
fast as the sorted-list algorithm, and asymptotically uses 
only one bit of memory per spin. 
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FIG. 8. The running times for the three algorithms for 
two-dimensional systems with R — 1.0 on a 266 MHz Pen- 
tium II with 128 MB of memory. Note that both the bits 
algorithm and the sorted-list algorithm have run times which 
grow approximately linearly (the logiV is not visible), and 
the brute force running time grows quadratically. Also notice 
that the largest bits simulation was 64 times larger than the 
largest sorted-list simulation. 

We have made CH — h source code implementing all 

three algorithms available at 

http://www.lassp.corncll.cdu/scthna/hystcrcsis/code/l . 
Each algorithm is contained in a separate, well- 
documented class. There are classes for detecting span- 
ning avalanches, measuring avalanche size distributions, 
and measuring correlation functions. There are also both 
Microsoft Windows and command line interfaces to the 
code. The command line interface should be portable 
to any computer with a C++ compiler and an imple- 
mentation of the Standard Template Library. We have 
also compiled executables for Windows 95/NT, Linux, 
and several other flavors of UNIX. The running times in 
Fig. H all come from the code compiled under Linux. 

Working on hysteresis, avalanches and noise with 
our model has been very rewarding. The simulations 
are beautiful and entertaining in themselves. Developing 
faster and more space efficient algorithms was amazingly 
satisfying: each new method not only eased our lives and 
our computer budgets, but also opened a whole new win- 
dow on the behavior of the model. It was also fun devel- 
oping new ways of measuring what was happening in the 
simulations: watching the spins flip and measuring the 
avalanche size distribution was only the beginning. We 
have written a few exercises which we hope will entertain 
and inform you as they did us. 

Suggested problems for further study 

1. Phase Transitions in the Shape of the Loop. Download 
our program from 

http: / /www. lassp.cornell.edu/scthna/hystercsis / code/ . 



Also download DynamicLattice and xmgr if you are using 
Unix. Run the program. Under Windows, just press OK 
in the opening dialog box to run a simulation with the 
default parameters. Under UNIX, type run at the >>> 
prompt. (More detailed instructions can be found on our 
Web site.) Try the other two algorithms. For smaller 
systems, brute force works acceptably, but for L = 500 
it is rather slow. 

You should see an animation of the avalanches as 
the external field is ramped upward and the spins flip. 
After the simulation ends, you should obtain a graph of 
M(H), the avalanche size distribution D(S,R), and the 
correlation function G(x,R). The M(H) curve shows 
the bottom half of the hysteresis loop: it should consist 
of many jumps of various sizes. The top half of the hys- 
teresis loop is pretty much the same shape —M{—H), 
but the details of the jumps are different. 

The avalanche size distribution D{S) measures the 
number of jumps as a function of size S. It should look 
like a fairly good power law, and be a straight line on a 
log-log plot. 

A log-log plot of the correlation function G(x, r) 
looks much less linear. At small distances, the corre- 
lation function decreases as a power law, but the power 
law behavior bends over after only a decade or so. This 
behavior is a symptom of R not being quite at the critical 
disorder R c . 

(a) According to our scaling theory P,[lO|], near the 
critical point D(S,R C ) ~ S~ T . What is your best esti- 
mate for the exponent f at the default value R = 1.0? 

(b) Do a simulation at a smaller value of disorder R, 
say 0.8. Does the behavior of G(x, R) remain a power 
law over larger values of x? 

(c) Do a simulation with D = 3, L = 50, and 
R = 2.5. What is f in three dimensions? (To obtain 
better data, do several runs and use the averaging op- 
tion.) Consider R = 2.1. Does the shape of the M(H) 
curve look qualitatively different? We believe that the 
power law distributions occur for R c ~ 2.16. At this 
value, the hysteresis loop first develops a macroscopic 
jump (the infinite avalanche). One obtains power laws 
and scaling at the phase transition in the shape of the 
hysteresis loop, between smooth loops and ones with an 
infinite avalanche. 

2. Time, Space, and Bits. We have a local, somewhat 
older supercomputer with four gigabytes of RAM. How 
large a system could we run in this memory using the 
three algorithms we discussed? Ignore all the memory 
requirements except those that scale linearly with the 
number of spins N; the rest is negligible. Assuming the 
time spent in bits continues to grow as N log N after the 
last data point in Fig. ||, how long will the largest possible 
bits simulation take to complete? 
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3. Programming. To do the following problems, you will 
need to download the source code for our program from 



http://www.lassp.corncll.cdu/scthna/hysteresis/code/l . 
You will also need a CH — h compiler which supports the 
standard template libraries. The standard template li- 
braries are an important advance for scientific program- 
ming. We find them incredibly useful, and look forward 
to their wide implementation now that the C++ stan- 
dard is in place. Links to compilers supporting the stan- 
dard template libraries can be found on our Web page. 

Run the simulation described in Problem f to test 
that the program works. Then try the following problems 
(many details on how to work with our source code can 
be found on our Web site). 

(a) Adding a Spin-Flip Action: Time Series for Large 
Avalanches. We have designed our code so that it is easy 
to add new types of measurements. One quantity which 
experimentalists measure is the change in magnetization 
with time. In our code, we record the time series of the 
whole run, with each avalanche represented as a single 
point. There is also interesting structure within each 
avalanche (see Fig. ||). Add a new class to the program 
which records the time series withinthe largest avalanche. 

(b) Implementing Brute Force. Implement the 
brute force algorithm. You can either replace the 
BruteForceHysteresisSimulation class or implement 
it from scratch. 

(c) Implementing Sorted 
Lists. Implement the sorted list algorithm by replacing 
the SortedListHysteresisSimulation class. 

(d) Implement the bits algorithm. 

More information on how 

to write your own BruteForceHysteresisSimulation, 
SortedListHysteresisSimulation, and 
BitsHysteresisSimulation classes are available on our 
Web site. 
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